2018
Você trabalha em uma empresa e a mesma não utiliza muitas técnicas para a previsão de demanda do produto “x”.
Após estudar previsão de demanda, você percebeu que as vendas deste produto possuem uma sazonalidade anual (o que é óbvio), porém a empresa não utiliza técnica de previsão de demanda.
Se a empresa possuísse uma previsão melhor, seria possível atuar de modo mais estratégico sobre o planejamento da produção (talvez antecipar a produção em meses de pico ao invés de realizar horas extra).
Você já pesquisou alguns softwares para previsão de demanda, porém ou são complicados de usar, ou possuem um preço proibitivo.
bibliotecas = c("forecast", "fpp2", "readxl")
install.packages(bibliotecas)
# Caso encontre algum erro, rode o install packages separadamente para cada biblioteca:
install.packages("forecast")
install.packages("fpp2")
install.packages("readxl")
library(forecast) library(fpp2) library(readxl)
# A série temporal "gold" foi carregada pela biblioteca forecast. autoplot(gold)
# Produção de Lâ na Austrália autoplot(woolyrnq)
# Produção de Gás na Austrália autoplot(gas)
# A função frequency determina a frequência da série. # Para dados sazonais, irá definir o período dominante da sazonalidade, # e para dados em ciclos, a duração média dos ciclos. frequency(gas)
## [1] 12
Observando a série temporal:
# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda Anti-Diabeticos")
Observando um gráfico sazonal
ggseasonplot(a10)
Observando um gráfico sazonal “polar”:
ggseasonplot(a10, polar = T)
Observando a Venda de Cerveja:
beer <- window(ausbeer, start = 1992)
autoplot(beer) + ylab("Venda de Cerveja")
Observando a Venda de Cerveja e sua Sazonalidade:
ggseasonplot(beer)
Mais uma forma de visualizar a demanda sazonal
ggsubseriesplot(beer)
library(readxl) dados_excel = readxl::read_xlsx(path = "dados/VendasCarros.xlsx", sheet = "Dados") head(dados_excel)
## # A tibble: 6 x 4 ## Ano Mês VendasCarros VendasCarrosAux ## <dbl> <dbl> <dbl> <dbl> ## 1 2018 1 2493 2486 ## 2 2018 2 2875 3003 ## 3 2018 3 3504 3454 ## 4 2018 4 3786 3665 ## 5 2018 5 4094 4353 ## 6 2018 6 4994 4842
ts_vendas = ts(data = dados_excel$VendasCarros, start = c(2018,1),frequency = 12) ts_vendas
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec ## 2018 2493 2875 3504 3786 4094 4994 4910 5575 5839 6202 6833 7382 ## 2019 2470 3059 3537 4003 4110 4516 5364 5854 5752 6089 7025 7801 ## 2020 2473 2884 3328 3689 4133 4932 5097 5397 6205 6374 6697 8024
forecast::autoplot(ts_vendas)
Previsão para a demanda \(y_{t+h}\), considerando dados disponiveis até \(y_{t}\)
Equação da Previsão: \[\hat{y}_{t+h|t} = \alpha y_{t} + \alpha(1-\alpha)y_{t-1} + \alpha(1-\alpha)^2y_{t-2} + ... \ para \ 0 \leq \alpha \leq 1\]
\[\hat{y}_{t+h|t} = \alpha y_{t} + \alpha(1-\alpha)y_{t-1} + \alpha(1-\alpha)^2y_{t-2} + ... \ para \ 0 \leq \alpha \leq 1\]
| Obs. | \(\alpha = 0.2\) |
|---|---|
| \(y_{t}\) | 0.2 |
| \(y_{t-1}\) | 0.16 |
| \(y_{-2}\) | 0.128 |
Previsão: \[\hat{y}_{t+h|t} = l_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)l_{t-1}\] O valor \(\alpha\) é obtido minimizando os erros ao rodar o modelo em um set de teste.
dados_petroleo = window(oil, start= 1996) previsao_petroleo = ses(dados_petroleo, h = 5) # Previsao para os próximos 5 anos summary(previsao_petroleo)
## ## Forecast method: Simple exponential smoothing ## ## Model Information: ## Simple exponential smoothing ## ## Call: ## ses(y = dados_petroleo, h = 5) ## ## Smoothing parameters: ## alpha = 0.8339 ## ## Initial states: ## l = 446.5868 ## ## sigma: 29.8282 ## ## AIC AICc BIC ## 178.1430 179.8573 180.8141 ## ## Error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set 6.401975 28.12234 22.2587 1.097574 4.610635 0.9256774 ## ACF1 ## Training set -0.03377748 ## ## Forecasts: ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 ## 2014 542.6806 504.4541 580.9070 484.2183 601.1429 ## 2015 542.6806 492.9073 592.4539 466.5589 618.8023 ## 2016 542.6806 483.5747 601.7864 452.2860 633.0752 ## 2017 542.6806 475.5269 609.8343 439.9778 645.3834 ## 2018 542.6806 468.3452 617.0159 428.9945 656.3667
autoplot(previsao_petroleo)
Previsão: \[\hat{y}_{t+h|t} = l_t + h b_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+b_{t-1})\] Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta)b_{t-1}\]
previsao_petroleo_holt = holt(dados_petroleo, h = 5) # Previsao para os próximos 5 anos autoplot(previsao_petroleo_holt)
summary(previsao_petroleo_holt)
## ## Forecast method: Holt's method ## ## Model Information: ## Holt's method ## ## Call: ## holt(y = dados_petroleo, h = 5) ## ## Smoothing parameters: ## alpha = 1e-04 ## beta = 1e-04 ## ## Initial states: ## l = 428.4833 ## b = 5.5771 ## ## sigma: 27.8684 ## ## AIC AICc BIC ## 177.2927 182.2927 181.7446 ## ## Error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set -0.2851768 24.57757 20.53231 -0.3285466 4.337459 0.8538816 ## ACF1 ## Training set 0.3429178 ## ## Forecasts: ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 ## 2014 534.4382 498.7234 570.1529 479.8172 589.0591 ## 2015 540.0148 504.3000 575.7295 485.3938 594.6357 ## 2016 545.5913 509.8766 581.3061 490.9704 600.2123 ## 2017 551.1679 515.4532 586.8827 496.5469 605.7889 ## 2018 556.7445 521.0298 592.4592 502.1235 611.3655
Previsão: \[\hat{y}_{t+h|t} = l_t + (\phi+\phi^2+...+\phi^h) b_t\]
Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+\phi b_{t-1})\]
Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) \phi b_{t-1}\] O parâmetro \(\phi\) é entre 0 e 1. Se o parâmetro for igual a 1, o crescimento será linear.
previsao_petroleo_holt = holt(dados_petroleo, h = 15, PI = F) previsao_petroleo_holt_damp = holt(dados_petroleo, h = 15, damped = T, PI = F) # Previsao para os próximos 5 anos autoplot(dados_petroleo) + autolayer(previsao_petroleo_holt, series="Holt Linear") + autolayer(previsao_petroleo_holt_damp, series="Holt com Damp")
Previsão: \[\hat{y}_{t+h|t} = l_t + hb_t + s_{t-m+h_m^+}\]
Nível: \[l_t = \alpha (y_t-s_{t-m}) + (1-\alpha)(l_{t-1}+b_{t-1})\]
Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]
Sazonalidade: \[b_t = \gamma(y_t - l_{t-1 - b_{t-1}}) + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a zero.
Previsão: \[\hat{y}_{t+h|t} = (l_t + hb_t) s_{t-m+h_m^+}\]
Nível: \[l_t = \alpha \frac{y_t}{s_{t-m}} + (1-\alpha)(l_{t-1}+b_{t-1})\]
Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]
Sazonalidade: \[b_t = \gamma \frac{y_t}{l_{t-1} + b_{t-1}} + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a 1.
Relembrando a série de anti-glicêmicos.
# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda AntiDiabeticos")
Relembrando a série de anti-glicêmicos.
# Produção de Gás na Austrália
previsao_anti_diab_aditivo = hw(a10, seasonal = "additive", PI= F)
previsao_anti_diab_multiplicativo = hw(a10, seasonal = "multiplicative", PI = F)
autoplot(a10) + ylab("Demanda de Remédios") +
autolayer(previsao_anti_diab_aditivo, series="HW Add.") +
autolayer(previsao_anti_diab_multiplicativo, series="HW Mult.")
previsao = naive(oil)
autoplot(oil, series="Dados") +
xlab("Ano") +
autolayer(fitted(previsao), series = "Previsao") +
ggtitle("Produção de Petróleo na Arábia Saudita")
autoplot(residuals(previsao))
checkresiduals(previsao)
## ## Ljung-Box test ## ## data: Residuals from Naive method ## Q* = 11.814, df = 9.8, p-value = 0.2824 ## ## Model df: 0. Total lags used: 9.8
treinamento = window(oil, end=2003)
teste = window(oil, start= 2004)
previsao = naive(treinamento,h = 10) # h = Número de períodos a prever
autoplot(previsao) +
ylab("Vendade Petróleo") +
autolayer(teste, series = "Dados de Teste")
accuracy(previsao, teste)
## ME RMSE MAE MPE MAPE MASE ## Training set 9.87358 52.56156 39.42504 2.506565 12.570647 1.0000000 ## Test set 21.60250 35.09832 29.97666 3.963914 5.777875 0.7603458 ## ACF1 Theil's U ## Training set 0.1801528 NA ## Test set 0.4029519 1.184862
library(Quandl)
ts_bovespa = Quandl("BCB/7", type = "ts")
autoplot(ts_bovespa)